Load packages

library(tidyverse)
library(readxl)
library(ggplot2)
library(zoo) # for function rollapply
library(scales) # for rescaling btw 0 and 1
library(plotly) # for interactive plot https://www.r-graph-gallery.com/163-interactive-area-chart-plotly.html
library(hrbrthemes) # for https://www.r-graph-gallery.com/163-interactive-area-chart-plotly.html 
library(htmlwidgets) # for https://www.r-graph-gallery.com/163-interactive-area-chart-plotly.html
library(viridis) # for https://www.r-graph-gallery.com/stacked-area-chart-plotly.html 
library(gganimate) #  for animated plot https://www.r-graph-gallery.com/287-smooth-animation-with-tweenr.html 
library(leaflet) # for map
library(magrittr) 
library(htmlwidgets)
library(rgdal)

Load data

#number of fishers
fishers<-read.csv(file.path("data/number_age_parttime_fulltime_fishers.csv"), sep=",")

#study municipalities list
municip<-read.csv(file.path("C:/Users/sigrid.engen/github/nor-prep/prep/administrative/komlist.csv"), sep=";")

#population in study municipalities
population<-read.csv(file.path("C:/Users/sigrid.engen/github/nor-prep/prep/administrative/data/population_size_study_municipalities_1986_2018.csv"),sep="," )

#spatial layer of municipalities in Norway
spat_municip<-readOGR(dsn="C:/Users/sigrid.engen/github/nor-prep/prep/administrative/raw/shapefiles/NO_AdminOmrader", "NO_AdminOmrader_pol")
## OGR data source with driver: ESRI Shapefile 
## Source: "C:\Users\sigrid.engen\github\nor-prep\prep\administrative\raw\shapefiles\NO_AdminOmrader", layer: "NO_AdminOmrader_pol"
## with 428 features
## It has 7 fields

Data cleaning - change municipal names (remove Norwegian letters and Sami names)

fishers_0<-left_join(fishers, municip, by = c("municip_numb" = "Komnum")) %>% 
  dplyr::select(-c(X, Kid, municip_name)) %>% 
  rename(municip_name = "Name")

Calculate 5 year moving average of number of fishers - main occupation

# summarise fishers across the age variable

fishers_1<- fishers_0 %>% 
  group_by(municip_numb,
           municip_name,
           year) %>%
  summarise(main_occup=sum(main_occupation)) %>% 
  ungroup()
## `summarise()` regrouping output by 'municip_numb', 'municip_name' (override with `.groups` argument)
# calculate moving average

### the rollapply function places the calculated value on the row = year three. 
### We want it to be placed year 6 - we want the moving average of the 5 previous years to
### compare this with the value the 6th year 

fishers_2 <- fishers_1 %>% 
  group_by(municip_numb) %>%
  mutate(five_yr_roll_aver = rollapply(main_occup, 5, mean, fill=NA)) %>% 
  dplyr::select(-main_occup) %>% 
  mutate(year = year + 3) %>%  #recoding year variable - have to remove count and add back in in next step so the order is correct
  filter(!year %in% c("2020", "2021", "2022")) %>% 
  ungroup()

# add values back to ensure that order of year and values are correct and calculate scores

fishers_3<- left_join(fishers_1, fishers_2, by = c("municip_numb", "municip_name", "year")) %>% 
  filter(!is.na(five_yr_roll_aver)) %>%     # remove NAs
  mutate(main_occup_score = main_occup/five_yr_roll_aver) %>%  #divide value the 6th year with average the 5 preceding years
  mutate(main_occup_score_1 = replace(main_occup_score, main_occup_score >= 1, 1)) # scores that are larger than 1 because the number of fishers are larger than the average the 5 preceeding years are given a value of 1

Plot scores - main occupation

ggplot(fishers_3)+
      geom_col(aes(x = year, y = main_occup_score_1), fill = "dodgerblue4")+
      facet_wrap(municip_name  ~ .) +
      labs(y = "Sustainability score", x = "") +
      ggtitle("")+
      theme(legend.position = "none", 
            axis.text.x = element_text(size = 5, angle = 90, hjust = 0.5, vjust = 0.5, face = "bold"),
            axis.text.y = element_text(size = 6, face = "bold", hjust = 0.5, vjust = 0.5, margin = margin(t = 0, r = 2, b = 0, l = 0)),
            plot.title = element_text(size=3, hjust=0.5, vjust=0.5),
            axis.ticks.x = element_line(size = 0.5),
            axis.ticks.length = unit(0.2, "cm"),
            text=element_text(family="sans"))

Calculate 5 year moving average of number of fishers - part-time occupation

fishers_1_p<- fishers_0 %>% 
  group_by(municip_numb,
           municip_name,
           year) %>%
  summarise(sum(part_time_occupation)) %>% 
  dplyr::rename(part_occup = 'sum(part_time_occupation)') %>% 
  ungroup()
## `summarise()` regrouping output by 'municip_numb', 'municip_name' (override with `.groups` argument)
fishers_2_p <- fishers_1_p %>% 
  group_by(municip_numb) %>%
  mutate(five_yr_roll_aver = rollapply(part_occup, 5, mean, fill=NA)) %>% 
  dplyr::select(-part_occup) %>% 
  mutate(year = year + 3) %>% 
  filter(!year %in% c("2020", "2021", "2022")) %>% 
  ungroup()

fishers_3_p<- left_join(fishers_1_p, fishers_2_p, by = c("municip_numb", "municip_name", "year")) %>% 
  filter(!is.na(five_yr_roll_aver)) %>% 
  mutate(part_time_occup_score = part_occup/five_yr_roll_aver) %>%
  mutate(part_time_occup_score_1 = replace(part_time_occup_score, part_time_occup_score >= 1 & part_time_occup_score < Inf, 1))

Plot scores - part time occupation

ggplot(fishers_3_p)+
      geom_col(aes(x = year, y = part_time_occup_score_1), fill = "dodgerblue4")+
      facet_wrap(municip_name  ~ .) +
      labs(y = "Number of fishers", x = "") +
      ggtitle("test")+
      theme(legend.position = "none", 
            axis.text.x = element_text(size = 5, angle = 90, hjust = 0.5, vjust = 0.5, face = "bold"),
            axis.text.y = element_text(size = 6, face = "bold", hjust = 0.5, vjust = 0.5, margin = margin(t = 0, r = 2, b = 0, l = 0)),
            plot.title = element_text(hjust=0.5, vjust=0.5),
            axis.ticks.x = element_line(size = 0.5),
            axis.ticks.length = unit(0.2, "cm"))
## Warning: Removed 5 rows containing missing values (position_stack).

Calculate 5 year moving average of number of fishers - all fishers

fishers_all<- fishers_0 %>% 
  group_by(municip_numb,
           municip_name,
           year) %>%
  summarise_at(c("main_occupation", "part_time_occupation"), sum) %>% 
  mutate(total_fishers = main_occupation+part_time_occupation) %>% 
  ungroup()

fishers_all_2 <- fishers_all %>% 
  group_by(municip_numb) %>%
  mutate(five_yr_roll_aver = rollapply(total_fishers, 5, mean, fill=NA)) %>% 
  select(-c(main_occupation:total_fishers)) %>% 
  mutate(year = year + 3) %>% 
  filter(!year %in% c("2020", "2021", "2022")) %>% 
  ungroup()

fishers_all_3<- left_join(fishers_all, fishers_all_2, by = c("municip_numb", "municip_name", "year")) %>% 
  filter(!is.na(five_yr_roll_aver)) %>% 
  mutate(occup_score = total_fishers/five_yr_roll_aver) %>%
  mutate(occup_score_1 = replace(occup_score, occup_score >= 1, 1))

Plot total number of fishers by municipality

ggplot(fishers_all_3)+
      geom_col(aes(x = year, y = total_fishers), fill = "dodgerblue4")+
      facet_wrap(municip_name  ~ .) +
      labs(y = "Number of fishers", x = "") +
      ggtitle("")+
      theme(legend.position = "none", 
            axis.text.x = element_text(size = 5, angle = 90, hjust = 0.5, vjust = 0.5, face = "bold"),
            axis.text.y = element_text(size = 6, face = "bold", hjust = 0.5, vjust = 0.5, margin = margin(t = 0, r = 2, b = 0, l = 0)),
            plot.title = element_text(hjust=0.5, vjust=0.5),
            axis.ticks.x = element_line(size = 0.5),
            axis.ticks.length = unit(0.2, "cm"))

### Plot total number of fishers by municipality (without Tromsø)

fishers_all_3_noTromso <- fishers_all_3 %>% 
  filter(!municip_name %in% c("Tromso"))

ggplot(fishers_all_3_noTromso)+
      geom_col(aes(x = year, y = total_fishers), fill = "dodgerblue4")+
      facet_wrap(municip_name  ~ .) +
      labs(y = "Number of fishers", x = "") +
      ggtitle("")+
      theme(legend.position = "none", 
            axis.text.x = element_text(size = 5, angle = 90, hjust = 0.5, vjust = 0.5, face = "bold"),
            axis.text.y = element_text(size = 6, face = "bold", hjust = 0.5, vjust = 0.5, margin = margin(t = 0, r = 2, b = 0, l = 0)),
            plot.title = element_text(hjust=0.5, vjust=0.5),
            axis.ticks.x = element_line(size = 0.5),
            axis.ticks.length = unit(0.2, "cm"))

Plot scores - all fishers

ggplot(fishers_all_3)+
      geom_col(aes(x = year, y = occup_score_1), fill = "dodgerblue4")+
      facet_wrap(municip_name  ~ .) +
      labs(y = "Sustainability score", x = "") +
      ggtitle("")+
      theme(legend.position = "none", 
            axis.text.x = element_text(size = 5, angle = 90, hjust = 0.5, vjust = 0.5, face = "bold"),
            axis.text.y = element_text(size = 6, face = "bold", hjust = 0.5, vjust = 0.5, margin = margin(t = 0, r = 2, b = 0, l = 0)),
            plot.title = element_text(hjust=0.5, vjust=0.5),
            axis.ticks.x = element_line(size = 0.5),
            axis.ticks.length = unit(0.2, "cm"))

### Plot distribution of all scores

ggplot(fishers_all_3, aes(occup_score_1)) +
  geom_bar(fill = "blue") +
  scale_x_binned()

### Plot distribution of scores by municipality

ggplot(fishers_all_3, aes(occup_score_1)) +
  facet_wrap(municip_name  ~ .) +
  geom_bar(fill = "blue") +
  scale_x_binned()+
  theme(axis.text.x = element_text(size = 5))

## Plot total number of fishers - interactive plot

fishers_northernNorw<- fishers %>% 
  group_by(year) %>%
  summarise_at(c("main_occupation", "part_time_occupation"), sum) %>% 
  mutate(total_fishers = main_occupation+part_time_occupation) %>% 
  ungroup()

# Usual area chart
p <- fishers_northernNorw %>%
  ggplot(aes(x=year, y=total_fishers)) +
  geom_area(fill="#69b3a2", alpha=0.5) +
  geom_line(color="#69b3a2") +
  ggtitle("Utvikling i antall fiskere i Nord-Norge 1983-2018,   Data: Fiskeridirektoratet") +
  ylab("Antall fiskere") +
  xlab("Ã…r") +
  theme_ipsum() +
  theme(plot.title = element_text(family = "arial", size=12, face= "italic"), 
        axis.title.x = element_text(family = "arial", size=10), 
        axis.title.y = element_text(family = "arial", size=10), 
        axis.line = element_line(size = 2, colour = "grey80"))

# Turn it interactive with ggplotly
p <- ggplotly(p)
p
# save the widget
#saveWidget(p, "figs/ggplotlyAreachart_Total_Fishers_NN.html", selfcontained = TRUE)

Plot total number of fishers - animated plot

fishers_northernNorw_long<- fishers_northernNorw %>% 
  select(-c(total_fishers)) %>%
  pivot_longer(cols = c("main_occupation", "part_time_occupation"),
                        names_to = "fishers", 
                        values_to = "n")

# Plot
p_2 <- fishers_northernNorw_long %>%
  ggplot(aes(x=year, y=n, group=fishers, color=fishers)) +
  geom_line() +
  geom_point() +
  scale_color_viridis(discrete = TRUE) +
  ggtitle("Utvikling i antall fiskere i Nord-Norge") +
  theme_ipsum() +
  ylab("Antall fiskere") +
  xlab("Ã…r")+
  transition_reveal(year) + 
  theme_set(theme_gray(base_size = 20, base_family = 'sans' )) +
  scale_color_manual(values=c('#008B8B','#000000'))
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
animate(p_2, duration = 5, fps = 20, width = 500, height = 500, renderer = gifski_renderer())

#anim_save("output.gif")

# Save at gif:
#anim_save(p_2, "figs/287-smooth-animation-with-tweenr.gif")
#getwd()

Plot fishers age distribution

fishers_age_NN<- fishers %>% 
  group_by(year, 
           age) %>%
  summarise_at(c("main_occupation", "part_time_occupation"), sum) %>% 
  mutate(total_fishers = main_occupation+part_time_occupation) %>% 
  ungroup()


#calculate age distribution in percentages
fishers_age_NN_percentages<- fishers_age_NN %>% 
  group_by(year) %>%
  mutate(percentage_main = main_occupation/sum(main_occupation)*100) %>% 
  mutate(percentage_part = part_time_occupation/sum(part_time_occupation)*100) %>%
  ungroup()

Plot fishers age distribution - main occupation

# Plot age distribution main occupation
p <- fishers_age_NN_percentages %>% 
  ggplot( aes(x=year, y=percentage_main, fill=age, text=age)) +
  geom_area( ) +
  scale_fill_viridis(discrete = TRUE) +
  theme(legend.position="none") +
  ggtitle("Aldersfordeling (%) blant fiskere med fiske som hovedyrke,   Data: Fiskeridirektoratet") +
  theme_ipsum() +
  theme(legend.position="right")+
  theme(plot.title = element_text(family = "arial", size=12, face= "italic"), 
        axis.title.x = element_text(family = "arial", size=10), 
        axis.title.y = element_text(family = "arial", size=10), 
        axis.line = element_line(size = 2, colour = "grey80"))

# Turn it interactive
p <- ggplotly(p, tooltip = c("text", "x", "fill"))
p
# save the widget
# library(htmlwidgets)
# saveWidget(p, file=paste0( getwd(), "/HtmlWidget/ggplotlyStackedareachart.html"))

Plot fishers age distribution - part time occupation

# Plot age distribution main occupation
p <- fishers_age_NN_percentages %>% 
  ggplot(aes(x=year, y=percentage_part, fill=age, text=age)) +
  geom_area( ) +
  scale_fill_viridis(discrete = TRUE) +
  theme(legend.position="none") +
  ggtitle("Aldersfordeling (%) blant fiskere med fiske som biyrke,   Data: Fiskeridirektoratet") +
  theme_ipsum() +
  theme(legend.position="right")+
  theme(plot.title = element_text(family = "arial", size=12, face= "italic"), 
        axis.title.x = element_text(family = "arial", size=10), 
        axis.title.y = element_text(family = "arial", size=10), 
        axis.line = element_line(size = 2, colour = "grey80"))

# Turn it interactive
p <- ggplotly(p, tooltip = c("text", "x", "fill"))
p
# save the widget
# library(htmlwidgets)
# saveWidget(p, file=paste0( getwd(), "/HtmlWidget/ggplotlyStackedareachart.html"))

Prepare data for mapping number of fishers

##################################
### Prepare attribute data #######
##################################

# Check correspondence btw the fishers file and study municip file in terms of municipal numbers
e<-anti_join(fishers, municip, by =c("municip_numb" = "Komnum")) # no municip numbers in fishers not in municip
f<-anti_join(municip, fishers, by =c("Komnum" = "municip_numb")) # noncoast is the only variable in municip not in fishers

# extract fishers data for 2018 for plotting on the map
fishers_northernNorw_2018<- fishers %>%
  filter(year =="2018") %>% 
  group_by(year, 
           municip_numb,
           municip_name) %>%
  summarise_at(c("main_occupation", "part_time_occupation"), sum) %>% 
  mutate(total_fishers = main_occupation+part_time_occupation) %>% 
  ungroup()

##################################
### Prepare spatial data #########
##################################

#Change municipal number of Harstad from 1901 to the correct 1903 in the spatial file over municipalities in NOrway. THe muncipality was merged with Bjarkjøy in 2013 and went from number 1901 to 1903. The polygon in the shapefile reflects this merger (Harstad includes Bjarkøy) but the number is wrong. 

spat_municip@data[spat_municip@data$KOMM == "1901", "KOMM"] <- "1903"

#extract subset of data (studymunicipalities from komlist.csv) from spatialpolygonsdataframe
muni_sub <- spat_municip[spat_municip$KOMM %in% c(1804L, 1805L, 1811L, 1812L, 1813L, 1815L, 1816L, 1818L, 1820L, 
          1822L, 1824L, 1827L, 1828L, 1832L, 1833L, 1834L, 1835L, 1836L, 
          1837L, 1838L, 1840L, 1841L, 1845L, 1848L, 1849L, 1850L, 1851L, 
          1852L, 1853L, 1854L, 1856L, 1857L, 1859L, 1860L, 1865L, 1866L, 
          1867L, 1868L, 1870L, 1871L, 1874L, 1902L, 1903L, 1911L, 1913L, 
          1917L, 1919L, 1920L, 1923L, 1924L, 1925L, 1926L, 1927L, 1928L, 
          1929L, 1931L, 1933L, 1936L, 1938L, 1939L, 1940L, 1941L, 1942L, 
          1943L, 2002L, 2003L, 2004L, 2012L, 2014L, 2015L, 2017L, 2018L, 
          2019L, 2020L, 2022L, 2023L, 2024L, 2025L, 2027L, 2028L, 2030L
          ), ]

#add number of fishers  to spatial layer of study municipalities
muni_fishers<- sp::merge(muni_sub, fishers_northernNorw_2018, by.x = "KOMM", by.y= "municip_numb")

#add names of municipalities from list  to spatial layer of study municipalities
muni_fishers_1<- sp::merge(muni_fishers, municip, by.x = "KOMM", by.y= "Komnum")

#transform projection of this spatial layer to fit leaflet and open street map
#https://stackoverflow.com/questions/34769845/how-to-plot-polygons-on-leaflet-using-r 
muni_fisher_transform <- spTransform(muni_fishers_1,  CRS("+ellps=WGS84 +proj=longlat +datum=WGS84 +no_defs"))


##################################
###Build map #####################
##################################

#https://www.theanalyticslab.nl/polygon-plotting-in-r/ 
#https://rstudio.github.io/leaflet/colors.html 
#https://www.r-graph-gallery.com/183-choropleth-map-with-leaflet.html


#Define cut points for the colorbins refering to number of fishers
cuts <- c(0.0, 1, 50, 100, 150, 200, 250, 300, 350, 400)

# Choose a color palette and assign it to the values
colorbins <- colorBin("Blues", domain =  muni_fisher_transform$total_fishers, bins = cuts)

# Prepare the text for tooltips:
mytext <- paste(
    "Kommune: ", muni_fisher_transform@data$Name,"<br/>", 
    "Antall fiskere: ", muni_fisher_transform@data$total_fishers) %>%
  lapply(htmltools::HTML)

# Show the number of fishers by municipality on a Leaflet map
map<-leaflet(muni_fisher_transform)  %>%
  setView(lng = 20, lat = 68, zoom = 5) %>%
  addTiles() %>% 
  addProviderTiles("Esri.WorldGrayCanvas") %>%
  addPolygons(stroke = TRUE, 
              color = "white", 
              weight="1", 
              smoothFactor = 0.3, 
              fillOpacity = 0.7, 
              fillColor = ~colorbins(muni_fisher_transform$total_fishers), 
              label = mytext,
              labelOptions = labelOptions(
                style = list("font-weight" = "normal", padding = "3px 8px"),
                textsize = "13px", 
                direction = "auto"
    )) 
  

# Add a legend
map_with_legend <- map %>% addLegend(pal = colorbins, 
                                      values = muni_fisher_transform$total_fishers,
                                      opacity = 0.7, title = "Number of fishers 2018", position = "bottomright")

map_with_legend
###########################
## SAVE OUTPUT HTML FILE ##
###########################

# Save output as HTML widget (or incorporate into Shiny / Flexdashboard)
#saveWidget(map_with_tooltip, file="Elderly residents in Utrecht.html")